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Abstract.  Driven  flow  of  a  non-equilibrium  non-conservative  (NENC)  system  with  a  mixture  of  immis¬ 
cible  particles  (A,  B  of  molecular  weight  Ma,Mb)  exhibits  self-organizing  patterns  (segregation,  phase- 
separation,  etc.)  in  steady-state.  The  flow  response  (v)  of  mass  flux  density  (j)  to  bias  (//),  v  =  dj/dH 
in  steady-state  is  found  to  be  sensitive  to  molecular  weight  ratio  (a  =  Mb/Ma )•  While  the  flux  density 
(j)  responds  linearly  to  bias  for  both  components  (A,  B)  at  a  =  1,  onset  of  eruptive  response  occurs  at 
extreme  bias  (H  — ►  1)  at  a  >  1  where  v  — >  oo  for  heavier  ( B )  and  v  — ►  — oo  for  lighter  (A)  constituents. 
Difference  in  molecular  weights  (Ma,Mb)  is  not  only  critical  to. eruptive  flow  but  also  in  controlling  the 
flow  response  prior  to  this  crossover. 

PACS.  61. 43. Bn  Structural  modeling:  serial-addition  models,  computer  simulation  83.10.Rs  Computer 
simulation  of  molecular  and  particle  dynamics 


Understanding  the  morphological  evolution,  transport, 
and  flow  [1-23]  of  particle  systems  have  attracted  enor¬ 
mous  attention  recently  due  to  their  response  to  external 
bias  particularly  in  granular  systems  [1-18].  Systems  mod¬ 
eled  by  coarse-grained  particles  have  diverse  variations  in 
spatial  and  temporal  scales,  i.e.,  nano-scale  materials  in 
laboratory  to  large  scale  applications  such  as  flow  of  wa¬ 
ter  and  sand  mixtures,  dissociation  of  methane  and  hydro¬ 
carbon  below  the  ocean  floor  [24-28],  and  mud  volcanoes. 
One  would  expect  interesting  responsive  phenomena  in  the 
flow  of  dissimilar  constituents  [29-31]  which  are  not  gener¬ 
ally  captured  by  traditional  hydrodynamic  methods.  Com¬ 
puter  simulations  with  particles  are  very  useful  in  under¬ 
standing  the  evolution  of  global  patterns  from  microscopic 
details.  A  range  of  computational  methods  are  employed 
in  such  studies  including  lattice  gas  [19,20],  molecular  dy¬ 
namics  [18,21,22],  and  Monte  Carlo  [23  methods.  Number 
of  particles  (constituents)  is  conserved  in  computational 
modeling  of  most  granular  systems  which  is  applicable  to 
a  range  of  systems  [1  18]  in  equilibrium.  In  many  driven 
systems  including  formation  and  dissociation  of  methane 
hydrate  and  mud  volcanoes,  the  number  of  particles  is  not 
conserved.  Non-equilibrium  steady-state  systems  [32,33] 
with  non-conserved  constituents  exhibit  interesting  flow 
response  and  structural  patterns  [29-31]. 


a  e-mail:  ras.pandeyOusm.edu 


Flow  properties  of  particulate  systems  are  frequently 
studied  by  a  number  of  lattice  gas  methods  [34  36].  It 
is,  howrever,  easier  to  incorporate  interactions  between 
constituent  particles  and  effective  medium  (empty  lattice 
sites,  see  below)  in  interacting  lattice  gas  [32,33]  than 
that  with  the  widely  used  Boltzmann  lattice  gas  [34,35] 
in  studying  the  flow.  Using  an  interacting  lattice  gas 
model  [29-31],  we  study  the  flow  response  of  an  immis¬ 
cible  mixture  with  dissimilar  components  A  and  B  driven 
by  a  pressure  bias  ( H )  against  sedimentation.  For  a  spe¬ 
cific  molecular  weight  ratio,  we  have  recently  observed  [31] 
an  eruptive  flow  response  at  the  extreme  value  of  the  bias. 
It  is  not  clear  how  such  a  quantity  as  the  difference  in 
molecular  weight  ratio  affect  the  response.  In  this  article 
we  show  that  the  difference  in  molecular  weights  of  dis¬ 
similar  constituents  is  not  only  key  to  eruptive  response 
but  also  in  orchestrating  the  flow  response  prior  to  its  di¬ 
vergence  at  the  extreme  bias. 

We  consider  a  mixture  of  twro  components  A  and  B 
(each  with  molecular  weight  Ma  and  Mb)  on  a  cubic  lat¬ 
tice  of  size  L3  with  L  =  30-200.  Almost  all  data  presented 
here  are  generated  on  1003  sample.  Different  sample  sizes 
are,  however,  used  to  check  for  the  finite  size  effects.  No  se¬ 
vere  finite  size  effect  is  observed  on  the  qualitative  behav¬ 
ior  unless  specified  otherwise.  The  molecular  weight  ratio 
o  =  Mb/M  a  is  varied  with  Ma  =0.1  and  a  =  1  —  10. 
Particles  A  and  B  are  randomly  distributed  initially  at 
about  50%  of  the  lattice  sites  with  one  particle  at  a  site. 
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A  nearest  neighbor  interaction  between  particles  (A,  B) 
and  empty  (pore)  sites  (O)  is  described  by  the  energy, 

^  =  ££J(*’n)  a) 

k  n 

where  index  k  runs  over  all  sites  occupied  by  particles 
and  n  over  all  nearest  neighbor  sites  of  /c.  The  interaction 
matrix  elements, 

J(A,  A)  =  J(B.  B)  =  -J(A.  B)  =  - J(B ,  A)  =  -e; 

J(A,0)  =  J(£,0)  -  -1.  (2) 


The  interaction  strength  or  miscibility  gap  6  =  1.  Molecu¬ 
lar  weights  of  constituents  affect  their  sedimentation.  The 
gravitational  potential  energy  Eg  of  a  particle  at  height  2 
(in  unit  of  g ,  the  acceleration  due  to  gravity)  is  given  by 

Eg  =  Ma/b  (3) 

The  sedimentation  probability  is  coupled  with  the  change 
in  gravitational  potential  energy  via  the  Boltzmann  dis¬ 
tribution  (see  below). 

A  source  of  particles  is  connected  to  the  bottom  plane 
(z  =  1).  This  causes  a  concentration  gradient  (see  below) 
which  exerts  a  driving  force  upward.  Additionally,  the 
hydrostatic  pressure  bias  (H)  pushes  fluid  constituents 
upward  (+z  direction)  against  the  gravitational  sedimen¬ 
tation  downward  (—z  direction).  The  bias  is  implemented 
in  selecting  one  of  the  nearest  neighbor  sites  along 
±yy  ±z  directions  with  probabilities, 


p  ~  p  —  p  —  p 

1  x  "  r—x  —  ry  1  -y  —  ^  i 


1  +  H  n  1 
Pz  =  - - - ,  P-z  =  ” 


H 


;  0  <  H  <  1. 


(4) 


Attempts  are  made  to  move  each  randomly  selected  parti¬ 
cle  ( A  and  B)  to  their  nearest  neighbor  sites  chosen  with 
the  bias  probability  H  with  the  Metropolis  algorithm.  Fol¬ 
lowing  procedure  is  used  to  implement  the  algorithm.  A 
particle  is  selected  randomly  say  at  site  i  and  one  of  its 
nearest  neighbor  site  j  is  selected  with  the  bias  probabil¬ 
ity  H.  If  the  site  j  is  already  occupied  by  another  particle, 
then  the  attempt  to  move  the  particle  fails.  If  the  site  j  is 
empty,  then  the  particle  is  moved  (from  the  site  i  to  site  j) 
with  the  Boltzmann  probability  e~AE/T  where  AE  is  the 
change  in  energy  E  =  Ei  4-  Eg  due  to  move  and  r  is  the 
temperature  in  unit  of  Boltzmann  constant  and  energy; 
r  —  1  is  used  in  this  study.  When  a  particle  at  the  bottom 
plane  (z  =  1)  moves  horizontally  (xy- plane)  or  upward  (to 
z  —  2  plane),  a  new  particle  (A  or  B)  from  the  source  is 
released  into  the  vacated  site  according  to  its  current  lat¬ 
tice  concentration  as  described  before  [29-31].  The  cubic 
box  is  open  along  vertical  boundaries,  i.e.,  a  particle  can 
drop  out  when  it  attempts  to  move  down  from  the  bot¬ 
tom  or  escape  from  the  top  (z  =  L).  Periodic  boundary 
conditions  are  implemented  along  the  transverse  (a:,  y)  di¬ 
rections.  An  attempt  to  move  each  particle  once  defines 
unit  Monte  Carlo  step  (MCS)  time. 

The  number  and  densities  of  particles  and  their  dis¬ 
tributions  change  as  the  simulation  proceeds.  Release  of 


particles  from  the  source  at  the  bottom  leads  to  a  net  flow 
along  the  longitudinal  (z)  direction.  The  competing  driv¬ 
ing  fields  due  to  pressure  bias  and  gravity  further  affect 
the  flow  and  pattern.  A  steady-state  is  reached  with  a  sta¬ 
ble  density  profile  and  a  constant  rate  of  mass  flux  in  the 
asymptotic  time  limit.  The  simulation  is  repeated  for  a 
number  of  independent  samples  to  obtain  a  reliable  esti¬ 
mate  of  physical  quantities.  These  quantities  include,  root 
mean  square  (rms)  displacement  of  tracers  (particles)  and 
that  of  their  center  of  mass,  density  profiles  (longitudinal 
and  transverse),  correlation  profiles  (i.e.,  average  number 
of  different  neighbors  of  each  particles),  flow  rate  and  flux 
rate  density,  etc.  as  a  function  of  the  bias  H  for  differ¬ 
ent  values  of  molecular  weight  ratios  a.  We  have  already 
studied  the  structural  response,  i.e.,  density  profiles  of  A 
and  B  and  their  steady-state  concentrations  [29].  In  this 
article,  we  constrain  to  flux  response. 

Figure  1  shows  a  typical  snap-shot  of  particles.  It  is 
clear  to  see  segregation  and  phase  separation  with  dense 
phase  toward  the  lower  part  and  low  density  toward  the 
top.  As  mentioned  above,  the  number  of  constituents  is 
not  conserved.  Therefore,  both  the  number  of  particles 
and  density  profiles  evolve  with  the  time  step  as  shown  in 
Figure  2  for  a  molecular  weight  ratio  a  =  10  at  high  values 
of  bias  H  =  0.7  1.0.  Initially  the  total  number  ( Na  )  of  A 
particles  in  the  lattice  decreases  while  the  number  (Nr)  of 
B  particles  increases  with  the  time  steps.  Both  Na  and  Nb 
attain  almost  constant  values  in  steady-state  equilibrium 
with  Nb  >  Na • 

The  equilibrium  values  of  Na  and  Nb  depend  non- 
monotonically  on  the  pressure  bias  [29].  The  steady-state 
density  (of  B  are  larger  than  that  of  A)  profiles  show 
non-linear  including  sigmoidal  variation  (see  below)  along 
the  longitudinal  direction  (from  bottom  to  top).  The  driv¬ 
ing  force  due  to  concentration  gradient  is  thus  non-linear 
along  the  longitudinal  direction  as  pointed  out  before  [29]. 
A  self-organizing  morphology  with  a  stable  structural  pat¬ 
tern  is  however  achieved  at  each  bias  and  molecular  weight 
ratio  even  though  there  is  a  constant  mobility  of  parti¬ 
cles  from  bottom  to  top.  At  high  pressure  bias  and  large 
molecular  weight  ratio  (a  =  10),  the  density  profile  is  over¬ 
whelmingly  dominated  by  the  heavier  (B)  particles  with 
a  relatively  uniform  distribution  along  the  transverse  di¬ 
rection  (Fig.  2). 

The  transverse  density  profiles  provide  insight  into 
the  complementary  structural  patterns.  For  example,  the 
spatial  complementary  oscillation  in  transverse  density 
profiles  of  A  and  B  reveals  a  phase  separation  or  de¬ 
segregation  [29]  at  low  bias  H  —  0.1  (Fig.  3).  Segrega¬ 
tion  (Fig.  3)  reduces  considerably  on  increasing  the  bias 
afid  almost  vanishes  at  high  values  of  the  bias  (Fig.  2) 
where  the  driving  bias  dominates  over  the  interaction  en¬ 
ergy  (E  =  Ei  - b  Eg).  The  difference  in  density  (higher 
with  higher  molecular  weight)  is  however  increases  with 
the  bias.  In  general,  increasing  the  pressure  bias  reduces 
the  non-linearity  [29]  in  the  longitudinal  density  profiles 
and  smears  out  the  contrast  (amplitude  of  oscillations)  in 
phase  separation  seen  from  the  transverse  density  profiles 
(Fig.  3).  These  structural  patterns  affect  the  flow  (follows). 
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Fig.  1.  Immiscible  A  (light,  Ma  =  0.1)  and  B  (dark,  Mb  =  0.3)  particles  with  pressure  bias  H  =  0.4  at  time  steps  t  =  1000 
on  a  253  lattice  with  e  =  1. 


Fig.  2.  Left:  number  of  particles  (Na,Nb)  versus  time  steps  at  bias  H  =  0.7  1.0  with  their  molecular  weights  Ma  =  0.1 
and  Mb  =  1.0.  A  (open,  black)  and  B  (filled)  symbols.  Right:  longitudinal  (z)  and  transverse  (y)  (inset)  density  profiles  in 
steady-state.  Sample  size  1003  with  128  independent  runs. 


Mass  Transfer 
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Fig.  3.  Density  profile  of  A 
(open)  and  B  (filled)  at  bias 
H  =  0.1  for  their  molecu¬ 
lar  weights  Ma  =  0.1  and 
Mb  =  0.1  —0.4.  Left:  longitu¬ 
dinal  (z)  and  right:  transverse 
( y )  density  profiles  in  steady- 
state.  Sample  size  1003  with 
128  independent  runs. 


Fig.  4.  Net  mass  transfer  along  longitudinal  direction 
of  A  and  B  (top  inset)  and  corresponding  flux  density 
(j)  versus  time  steps  at  bias  H  —  1.0  for  Ma  =  0.1  and 
Mb  =  0.1, 0.3,  0.4,  0.8, 1.0.  As  many  as  128  independent 
samples  are  used  on  1003  lattice. 


Obviously,  there  is  a  net  flow  of  A  and  B  from  bottom 
to  top.  Evolution  of  the  mass  flux  with  the  time  step  is 
presented  in  Figure  4.  The  mass  transfer  Q  of  each  compo¬ 
nent  increases  linearly  with  the  time  steps  in  the  asymp¬ 
totic  long  time  regime  which  indicates  that  the  system 
has  reached  a  steady-state.  Corresponding  current  (flux) 
density  j  can  be  evaluated  from 

1  dQ 
3  Lx  X  Ly  dt  ' 


Evolution  of  j  shows  howr  it  reaches  a  constant  value  for 
each  molecular  weight  ratio.  The  constant  value  of  the  flux 
density  implies  that  the  dynamic  system  has  reached  the 
steady-state  wdiich  is  also  consistent  with  a  stable  mor¬ 
phology.  Thus,  both  structure  and  flow  become  stable  in 
the  asymptotic  time  regime.  It  is  important  to  point  out 
that  the  unit  of  time  is  arbitrary  similar  to  energy  and 
molecular  weight  which  is  a  common  practice  in  such  a 
coarse-grained  modeling  with  phenomenological  interac¬ 
tions.  However,  the  time  step  appears  to  scale  linearly 


R.B.  Pandey  and  J.F.  Gettrust:  Eruption  in  a  driven  flow  of  self-organizing  immiscible  system 


87 


Fig.  5.  Steady-state  flux  density  (j)  of  A  and  B 
versus  bias  H  for  Ma  =0.1  and  Mb  =  0.1  1.0.  As 
many  as  128  independent  samples  are  used  on  1003 
lattice. 


with  a  realistic  time  frame  born  out  by  the  linear  mass 
flow  expected  from  the  solution  of  a  simple  driven  dif¬ 
fusive  system  [37],  Our  goal  is  to  understand  trends  in 
the  flow  response  with  respect  to  variations  in  parameters 
such  as  pressure  bias  and  molecular  weight. 

How  does  the  steady-state  flux  density  (Ja/b)  re¬ 
spond  to  pressure  bias  for  mixtures  of  different  molec¬ 
ular  weight  ratios?  Figure  5  shows  the  variation  of  the 
flux  density  with  the  pressure  bias  for  various  mixtures, 
a  =  Mb/Ma  =  1-10.  For  a  symmetric  molecular  weight 
ratio  (a  =  1),  we  see  a  linear  response  throughout  the 
pressure  bias  regime.  On  increasing  the  molecular  weight 
ratio,  the  flux  density  response  of  A  and  B  begin  to  dif¬ 
fer.  At  low  bias  the  flux  density  response  of  both  com¬ 
ponent  is  linear  (which  is  better  seen  in  the  normal  plot, 
the  inset  figure  than  on  the  log-log  plot).  At  large  values 
of  bias,  the  response  becomes  non-linear.  The  response  of 
higher  molecular  weight  (Mb  —  0.2  1.0)  component  di¬ 
verges  while  that  of  lower  weight  component  (A)  declines 
dramatically  at  extreme  values  of  the  bias.  Let  us  define 
a  volatility  index, 

v  =  dj/dH.  (6) 

We  see  that  v  — ►  oo  for  B  and  v  — ►  — oc  for  A  as  H  — >  1. 

What  happens  when  a  difference  in  the  molecular 
weight  of  particles  A  and  B  sets  in  that  causes  an  erup¬ 
tive  flow  response  of  the  heavier  (B)  component?  We  have 
already  seen  such  an  eruptive  response  for  a  fixed  value  of 
q  =  3  [31].  The  variation  of  the  molecular  weight  ratio  (a) 
as  presented  here,  however,  shows  that  (?)  the  difference 
in  molecular  weight  is  the  root  cause  of  eruptive  response 
(which  vanishes  for  a  —  1),  (n)  approach  to  eruptive  flow 
differ  considerably  with  the  difference  in  molecular  weight 


of  the  constituents  larger  the  difference,  the  faster  is  the 
response  (v)  of  the  heavier  component  prior  to  eruption, 
and  (iii)  the  nature  of  eruption  remains  independent  of 
the  molecular  weight  ratio  (a  >  1)  as  H  — ►  1.  Asym¬ 
metry  in  the  molecular  weights  of  the  constituents  seems 
to  introduce  correlation  between  the  source  and  the  dy¬ 
namic  lattice  system.  Particularly,  the  correlation  among 
particles  with  the  higher  molecular  weight  increases  on  in¬ 
creasing  the  bias  leading  to  eruption  in  their  flow  while  an 
opposite  trend  occurs  concurrently  in  flow  of  the  compo¬ 
nent  with  the  lower  molecular  weight. 
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